I. Table1

Demographic characteristics of our sample depending on the group.

dat_wide$age <- as.numeric(as.character(dat_wide$age))
dat_wide$sex <- factor(dat_wide$sexe)

table1::table1(~ age + sex | group, data=dat_wide, overall = "Overall",
              render.continuous=c(.="Mean / Median<br>SD<br>[min, max]"))
Control
(N=22)		 Experimental
(N=20)		 Overall
(N=42)
age
Mean / Median
SD
[min, max]		 63.1 / 62.0
3.64
[58.3, 69.5]		 63.8 / 64.2
3.33
[57.8, 69.5]		 63.4 / 63.3
3.47
[57.8, 69.5]
sex
F 12 (54.5%) 11 (55.0%) 23 (54.8%)
M 10 (45.5%) 9 (45.0%) 19 (45.2%)

Effect sizes of the group differences

smd_age <- with(dat_wide, lsr::cohensD(age ~ group))
smd_age
## [1] 0.2111493
chi <- with(dat_wide, chisq.test(table(sexe, group)))
confintr::cramersv(chi)
## [1] 0

II. Distributions

Outcome distribution at baseline

Inhibition

dat_inhib_dist <- subset(dat_long, grepl("Inhibition", dat_long$GEN_outcome) & time != "t1")

ggplot(dat_inhib_dist, aes(x = value)) + 
  geom_histogram(aes(y=..density..), fill = "grey", alpha = 0.8) + 
  geom_density(alpha = 0.3, aes(fill = group)) + 
  facet_wrap(. ~ outcome, scales = "free") +
  theme_bw() + 
  theme(legend.position = "bottom")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Working Memory

dat_memory_dist <- subset(dat_long, grepl("memory", dat_long$GEN_outcome) & time != "t1")

ggplot(dat_memory_dist, aes(x = value)) + 
  geom_histogram(aes(y=..density..), fill = "grey", alpha = 0.8) + 
  geom_density(alpha = 0.3, aes(fill = group)) + 
  facet_wrap(. ~ outcome, scales = "free") +
  theme_bw() + 
  theme(legend.position = "bottom")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Overal EF

dat_EF_dist <- subset(dat_long, grepl("Executive", dat_long$GEN_outcome) & time != "t1")

ggplot(dat_EF_dist, aes(x = value)) + 
  geom_histogram(aes(y=..density..), fill = "grey", alpha = 0.8) + 
  geom_density(alpha = 0.3, aes(fill = group)) + 
  facet_wrap(. ~ GEN_outcome, scales = "free") +
  theme_bw() + 
  theme(legend.position = "bottom")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Other

dat_othr_dist <- subset(dat_long, !grepl("memory", dat_long$GEN_outcome) & !grepl("Inhibition", dat_long$GEN_outcome) & !grepl("Executive", dat_long$GEN_outcome) & time != "t1")

ggplot(dat_othr_dist, aes(x = value)) + 
  geom_histogram(aes(y=..density..), fill = "grey", alpha = 0.8) + 
  geom_density(alpha = 0.3, aes(fill = group)) + 
  facet_wrap(. ~ outcome, scales = "free") +
  theme_bw() + 
  theme(legend.position = "bottom")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

III. Preliminary analyses

BeeBot mastering

dat_eval = subset(dat_wide, s4_bin == 1)
paste0(round(nrow(dat_eval)/sum(!is.na(dat_wide$s4_bin))*100, 3), "% of the sample succeed at a beebot mastering test after the last training session, with a mean score of ", round(mean(dat_eval$s4_0_4), 3), "/4")
## [1] "89.474% of the sample succeed at a beebot mastering test after the last training session, with a mean score of 3.824/4"

Correlation between tests

ggstatsplot::ggcorrmat(dat_wide1[,which(colnames(dat_wide1) == "Denomination (NEPSY-II)"):ncol(dat_wide1)],
          p.adjust.method = "none",
          ggtheme = ggplot2::theme_bw(),
          colors = c("#0016FF", "white", "#FF0000"),
          sig.level = 0.99)

# ggstatsplot::ggcorrmat(dat_wide2[,which(colnames(dat_wide2) == "Planification (KABC-II)"):ncol(dat_wide2)],
#           p.adjust.method = "none",
#           ggtheme = ggplot2::theme_bw(),
#           colors = c("#0016FF", "white", "#FF0000"),
#           sig.level = 0.05)

IV. Main analysis

Table

res_main = res_main %>% mutate(across(where(is.numeric), round, digits=3))
DT::datatable(res_main, 
              rownames = FALSE,
              extensions = "Buttons",
              options = list(  # options
                scrollX = TRUE,
                dom = c('tB'), 
                buttons = c('copy', 'csv', 'excel', 'pdf', 'print'),
                scrollY = "600px", 
                pageLength = 30,
                order = list(4, 'asc'),
                columnDefs = list(
                  list(width = '70px',
                       targets = "_all"),
                  list(className = 'dt-center', 
                       targets = "_all"))))

Forest plot

tab.main <- data.frame(
  'General outcome' = res_main$outcome,
  TOST = res_main$TOST)

metaviz::viz_forest(x = res_main[, c("d", "d_se")], 
           variant = "classic",
           col = "Greys", xlab = "SMD", annotate_CI = TRUE, 
           study_table = tab.main,
           type = "study_only",
           text_size = 2.5,
           x_limit = c(-1.5, 1.5),
           x_breaks = seq(-3, 3, 1))

Box plots

Inhibition

dat_inhib <- subset(dat_plot, grepl("Inhibition", dat_plot$GEN_outcome))

ggplot(dat_inhib, aes(x = time, y = value, color = group, fill = group)) +
  # geom_line() +
  geom_boxplot(alpha = 0.2) +
  facet_grid(. ~ GEN_outcome,  scales = 'free') + 
theme_bw() + 
  theme(legend.position="bottom")

Working memory

dat_memory <- subset(dat_plot, grepl("memory", dat_plot$GEN_outcome))

ggplot(dat_memory, aes(x = time, y = value, color = group, fill = group)) +
  geom_boxplot(alpha = 0.2) +
  facet_grid(. ~ GEN_outcome,  scales = 'free') + 
theme_bw()  + 
  theme(legend.position="bottom")

Overall EF

dat_EF <- subset(dat_plot, grepl("Executive", dat_plot$GEN_outcome))

ggplot(dat_EF, aes(x = time, y = value, color = group, fill = group)) +
  geom_boxplot(alpha = 0.2) +
  facet_grid(. ~ GEN_outcome,  scales = 'free') + 
theme_bw()  + 
  theme(legend.position="bottom")

Others

dat_othr <- subset(dat_plot, !grepl("memory", dat_plot$GEN_outcome) & !grepl("Inhibition", dat_plot$GEN_outcome) & !grepl("Executive", dat_plot$GEN_outcome))

ggplot(dat_othr, aes(x = time, y = value, color = group, fill = group)) +
  geom_boxplot(alpha = 0.2) +
  facet_grid(. ~ GEN_outcome,  scales = 'free') + 
  theme_bw() + 
  theme(legend.position="bottom")

V. Additional analyses

1. No demographic covariates

## Warning in qnorm(1 - cer): NaNs produced

Table

res_S2 = res_S2 %>% mutate(across(where(is.numeric), round, digits=3))
DT::datatable(res_S2, 
              rownames = FALSE,
              extensions = "Buttons",
              options = list(  # options
                scrollX = TRUE,
                dom = c('tB'), 
                buttons = c('copy', 'csv', 'excel', 'pdf', 'print'),
                scrollY = "600px", 
                pageLength = 30,
                order = list(4, 'asc'),
                columnDefs = list(
                  list(width = '70px',
                       targets = "_all"),
                  list(className = 'dt-center', 
                       targets = "_all"))))

Forest plot

tab.S2 <- data.frame(
  'General outcome' = res_S2$outcome,
  TOST = res_S2$TOST)

metaviz::viz_forest(x = res_S2[, c("d", "d_se")], 
           variant = "classic",
           col = "Greys", xlab = "SMD", annotate_CI = TRUE, 
           study_table = tab.S2,
           type = "study_only",
           text_size = 2.5,
           x_limit = c(-1.5, 1.5),
           x_breaks = seq(-3, 3, 1))

2. All outcomes

## Warning in qnorm(1 - cer): NaNs produced

## Warning in qnorm(1 - cer): NaNs produced

## Warning in qnorm(1 - cer): NaNs produced

Table

res_S1 = res_S1 %>% mutate(across(where(is.numeric), round, digits=3))

DT::datatable(res_S1, 
              rownames = FALSE,
              extensions = "Buttons",
              options = list(  # options
                scrollX = TRUE,
                dom = c('tB'), 
                buttons = c('copy', 'csv', 'excel', 'pdf', 'print'),
                scrollY = "600px", 
                pageLength = 30,
                order = list(4, 'asc'),
                autoWidth = TRUE,
                columnDefs = list(
                  list(width = '70px',
                       targets = "_all"),
                  list(className = 'dt-center', 
                       targets = "_all"))))

Forest plot

tab.supp <- data.frame(
  Outcome = res_S1$outcome,
  TOST = res_S1$TOST)

metaviz::viz_forest(x = res_S1[, c("d", "d_se")], 
           variant = "classic",
           col = "Greys", xlab = "SMD", annotate_CI = TRUE, 
           study_table = tab.supp,
           type = "study_only",
           text_size = 2.2,
           x_limit = c(-1.5, 1.5),
           # N = tab.prim$N,
           x_breaks = seq(-3, 3, 1))

3. Beebot mastering

dat_semiwide$z_change = dat_semiwide$z1 - dat_semiwide$z0
dat_s3 = subset(dat_semiwide, !is.na(s4_0_4))
dat_s3$s4_bin = factor(dat_s3$s4_bin)

summary(lmerTest::lmer(z_change ~ ordered(s4_0_4) + (1|participant) + (1|outcome), data = dat_s3))
## boundary (singular) fit: see ?isSingular
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: z_change ~ ordered(s4_0_4) + (1 | participant) + (1 | outcome)
##    Data: dat_s3
## 
## REML criterion at convergence: 640.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2475 -0.4357 -0.0847  0.4683  3.3823 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  participant (Intercept) 0.03073  0.1753  
##  outcome     (Intercept) 0.00000  0.0000  
##  Residual                0.93392  0.9664  
## Number of obs: 228, groups:  participant, 19; outcome, 12
## 
## Fixed effects:
##                    Estimate Std. Error        df t value Pr(>|t|)
## (Intercept)        0.064602   0.131962 15.000001   0.490    0.632
## ordered(s4_0_4).L  0.003478   0.245454 15.000001   0.014    0.989
## ordered(s4_0_4).Q -0.061161   0.263923 15.000001  -0.232    0.820
## ordered(s4_0_4).C  0.072150   0.281182 15.000001   0.257    0.801
## 
## Correlation of Fixed Effects:
##             (Intr) o(4_0_4).L o(4_0_4).Q
## or(4_0_4).L -0.618                      
## or(4_0_4).Q -0.169 -0.431               
## or(4_0_4).C  0.093 -0.102     -0.398    
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
res_s3_bin = lmerTest::lmer(z_change ~ s4_bin + (1|participant) + (1|outcome), data = dat_s3)
## boundary (singular) fit: see ?isSingular
summary(res_s3_bin)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: z_change ~ s4_bin + (1 | participant) + (1 | outcome)
##    Data: dat_s3
## 
## REML criterion at convergence: 637.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2931 -0.4380 -0.0677  0.4747  3.4470 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  participant (Intercept) 0.01485  0.1219  
##  outcome     (Intercept) 0.00000  0.0000  
##  Residual                0.93392  0.9664  
## Number of obs: 228, groups:  participant, 19; outcome, 12
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)
## (Intercept)   0.2217     0.2153 17.0000    1.03    0.317
## s4_bin1      -0.1866     0.2276 17.0000   -0.82    0.424
## 
## Correlation of Fixed Effects:
##         (Intr)
## s4_bin1 -0.946
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
emmeans(res_s3_bin, ~s4_bin)
##  s4_bin emmean     SE    df lower.CL upper.CL
##  0      0.2217 0.2153 16.79   -0.233    0.676
##  1      0.0351 0.0738  8.85   -0.132    0.203
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95
ggplot(dat_s3, aes(x = s4_0_4, y = z_change)) + 
  geom_jitter(width=0.2, size = 2, alpha = 0.4) + 
  theme_bw()

ggplot(dat_s3, aes(x = factor(s4_bin), y = z_change)) + 
  geom_jitter(width=0.1, size = 2, alpha = 0.4) + 
  theme_bw()

# dat_agg <-
#   dat_long %>%
#   # pivot_wider(names_from = "time", values_from = "value") %>%
#   dplyr::group_by(outcome) %>%
#   dplyr::mutate(z = scale(value)) %>%
#   dplyr::ungroup() %>%
#   dplyr::group_by(time, ID, group) %>%
#   summarise(value_z = sum(z))
# 
# ID <- 
#   dat_agg %>%
#   group_by(ID) %>%
#   arrange(time) %>%
#   mutate(diff = value_z - lag(value_z, default = first(value_z))) %>%
#   filter(time == "t1") %>%
#   mutate(Improvement = if_else(diff > 0, "Improvement", "Deterioration"))%>%
#   dplyr::select(ID, Improvement) %>%
#   filter(!is.na(Improvement))
#  
# dat_line <- merge(dat_agg, ID, by = "ID")
# ggplot(dat_line, aes(x = time, y = value_z, color = group, group = ID)) +
#   geom_point(size = 2.5, alpha = 0.4) + 
#   geom_line(size = 1, linetype  = "dotted", alpha = 0.8) +
#   theme_bw() + 
#   facet_grid(group ~  Improvement) +
#   theme(legend.position = "none")